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The grant period covered in this report was devoted to the study of radiative transfer 
approximations in collaboration with Dr. Michael King of NASA/Goddard Space Flight Center. 
A manuscript describing the study has been submitted for publication to the Journal of the 
Atmospheric Sciences. A copy is being attached for your records. 

In addition the prior study on surface longwave fluxes derived using a hybrid method 
(reported earlier) has been re-submitted after revisions to the Journal of Climate. The revised 
manuscript is also attached for your records. 

As part of the ongoing research project we have now acquired a continuous record of 
ISCCP C-l data extending from October 1, 1985 - April 30, 1989. This extensive cloud data set 
will be used in future studies including efforts in collaboration with personnel in the Radiation 
and Climate Branch of NASA/Goddard Space Flight Center. 

Some ongoing research results will be presented at the AGU Spring meeting in Montreal, 
May 11-16. The titles of two presentations to be made at the meeting are: 

1 . The application of similarity relations to the computation of spectrally integrated 
solar absorption by water clouds by Harshvardhan and Michael D. King 
(NASA/Goddard). 

2. A spectral investigation of the atmospheric shortwave cloud radiative forcing by 
William L. Ridgway (ARC; NASA/Goddard) and Harshvardhan. 


CLOUD RADIATIVE FORCING AND THE SURFACE 
LONGWAVE RADIATION 


Hui Zhi and Harshvardhan 


Department of Earth and Atmospheric Sciences 
Purdue University 
West Lafayette, IN 47907 


Abstract 


Global maps of the monthly mean net upward longwave radiation flux at the ocean 
surface have been obtained for April, July, October 1985 and January 1986. These maps were 
produced by blending information obtained from a combination of general circulation model 
cloud radiative forcing fields, the top-of-the-atmosphere cloud radiative forcing from ERBE and 
TOVS profiles and sea surface temperature on ISCCP Cl tapes. The fields are compatible with 
known meteorological regimes of atmospheric water vapor content and cloudiness. There is a 
vast area of high net upward longwave radiation flux (> 80 Wm' 2 ) in the eastern Pacific Ocean 
throughout most of the year. Areas of low net upward longwave radiation flux (<40 Wm -2 ) are 
the tropical convective regions and extra tropical regions that tend to have persistent low cloud 
cover. The technique used in this study relies on GCM simulations and so is subject to some of 
the uncertainties associated with the model. However, all input information regarding 
temperature, moisture and cloud cover is from satellite data having near global coverage. This 
feature of the procedure alone warrants its consideration for further use in compiling global maps 
of the net longwave radiation at the surface, at least over the oceans. 
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1. Introduction 

The surface radiation budget, i.e., the net solar radiation absorbed minus the net 
longwave radiation emitted, and its spatial and temporal variations are key parameters in climate 
and weather studies. This budget plays a major role in determining radiative heating, as well as 
sensible and latent heat fluxes over ocean and land surfaces. As a result, the net radiative flux 
constitutes an important boundary forcing for the general ocean circulation and a crucial 
parameter for determining meridional oceanic heat transport, ocean-atmosphere interaction and 
land-atmosphere interaction. Moreover, it is a useful parameter when addressing issues related 
to climate change due to CO 2 and other trace gases, and in the validation of radiation schemes 
used in climate models. Therefore, it is understandable for the atmospheric and oceanic 
communities to need reliable estimates of the surface radiation budget (WCP-92 1984). 

Direct high-quality radiation measurements at the surface are difficult to make, 
particularly over the oceans which cover more than 60% of the Earth's surface. Actually there 
are very few surface stations measuring the radiation budget routinely and reliably because of the 
requirement of careful instrument calibration and temperature correction for the radiation, 
especially longwave, measurement. In addition, because of operating costs, it is not feasible to 
maintain a network of surface stations over the oceans. Although an attempt has been made to 
use ships to observe some meteorological parameters such as sea surface temperature, air 
temperature, specific humidity near the surface and even fractional cloud coverage, there has not 
been much progress in the measurement of surface radiation. Few ships measure radiation 
quantities because of special needs that require a dedicated facility for this purpose. Moreover, 
regular ship observations are limited along commercial shipping lanes, and vast geographic gaps 
still exist, especially in the southern hemisphere. Consequently, the empirical formulas used to 
derive budgets have been validated only over limited regions, and when applied globally, large 
errors are inevitable. Therefore, direct measurement of shortwave and longwave radiation fluxes 
at the surface globally has not been possible. 

Since there are difficulties in obtaining radiative data from surface stations routinely and 
reliably, it has been realized that space based observations are the only means to have global 
coverage. However, because of the intervening atmosphere, the surface radiation budget is 
difficult to measure from satellites, whereas the top-of-the-atmosphere (TOA) radiation balance 
can be measured directly. Over the past decades, considerable effort has been expended in the 
global measurement of the TOA radiation budget. Attempts at inferring the surface radiation 
budget from space based measurements have only begun recently. 

There has been some success in obtaining the solar radiation budget at the surface (e.g. 
Raschke and Preuss 1979; Tarpley 1979; Gautier et al. 1980; Pinker and Ewing 1985; Justus et 
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al. 1986). The progress in surface longwave radiation budget measurements from space, 
however, has been much slower. Currently, some techniques are available to estimate the 
downward longwave component of the surface flux (Darnell et al. 1983, 1986; Chou 1985; 
Schmetz et al. 1986; Frouin et al. 1988; Ardanuy et al. 1989; Wu and Cheng 1989; Breon et al. 
1991). The net longwave component at the surface can be estimated by the difference between 
the upwelling and downwelling fluxes. The upward component is determined directly from sea 
surface temperature since the oceanic surface emits essentially as a blackbody. The downward 
flux, however, is more difficult to obtain since it depends on many meteorological parameters 
such as atmospheric moisture, temperature and cloud cover. Because of the uncertainty of the 
measurement of these meteorological parameters, at present there is a need for improvement in 
the estimation of net longwave radiation fluxes at the surface. 

Two types of methods have been used to estimate the downward longwave flux at the 
surface: statistical and physical. As the name implies, statistical methods Tely on correlations 
between fluxes and observed meteorological parameters. The physical techniques are based on 
modeling radiative processes occurring in the atmosphere (clear and cloudy atmosphere). The 
downward flux is computed from radiative transfer models which utilize parameters obtained 
from satellite radiance data. These parameters include temperature and water vapor mixing ratio 
profiles, fraction of cloud coverage and cloud emittance. However, all physical methods 
currently under consideration have to make certain assumptions regarding both the presence of 
clouds and their vertical extent. Recent examples of these attempts are Chou (1985), Schmetz et 
al. (1986), Darnell et al. (1986), Gupta (1989), and Wu and Cheng (1989). The treatment of 
longwave radiation transmittance in the presence of clouds becomes more complex since 
knowledge of cloud top and base heights and emittances are required. For this reason, it is 
important to determine the vertical profile of cloudiness as well as the horizontal distribution of 
clouds and associated emittances. Unfortunately, determination of the vertical profile of 
cloudiness from space based measurements is difficult since overlap of cloudy layers is common 
in the real atmosphere. Therefore, because of the uncertainties in assumed cloudiness, all these 
methods often give unreliable results. 

The method used here to obtain monthly mean quantities avoids the explicit computation 
of cloud fraction and the location of cloud base in estimating the downward longwave radiation 
globally (Harshvardhan et al. 1990). An advantage of this technique is that no independent 
knowledge or assumptions regarding cloud cover for a particular month are required. The only 
infonnation required is a relationship between the cloud radiative forcing (CRF) at the top of the 
atmosphere and that at the surface, which is obtained from a general circulation model (GCM) 
simulation. 
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2. Method 

The cloud radiative forcing (CRF) has been defined as the difference between the 
radiative flux as measured or computed and the clear sky flux (Charlock and Ramanathan 1985; 
Ramanathan 1987). For example, the longwave cloud radiative forcing at the surface is 

Surface LWCRF = Surface LW Flux - Clear Sky Surface LW Flux. (1) 

Previous studies have shown that there is no correlation between the spectrally integrated 
outgoing longwave radiation and the net longwave at the surface (Ramanathan 1986; Weare 
1989). But Harshvardhan et al. (1990) have recently come to the conclusion that there is a 
relationship between the longwave CRF at the top of the atmosphere and the surface in model 
simulations. This relationship has been explained on the grounds that clouds of a certain type 
tend to form preferentially over certain geographic areas. The technique proposed by 
Harshvardhan et al. (1990) to estimate surface longwave fluxes starts with the GCM simulated 
climatological ratio of the CRF at the top and at the surface along with the longwave CRF at the 
top of the atmosphere obtained from the Earth Radiation Budget Experiment (ERBE; 
Ramanathan 1987) to compute the surface longwave CRF for the particular month. In this way, 
no independent knowledge or assumptions about cloud cover are involved, thus avoiding the 
most uncertain step in other methods of estimating the longwave radiation budget at the surface. 
The next step in the procedure is to obtain an estimate of the clear sky downward longwave flux 
and upward emission at the surface. 

As mentioned before, the downwelling longwave radiation at the surface can not be 
measured directly from space, but profiles of temperature and water vapor mixing ratio in clear 
columns are routinely obtained from inversions of measured radiances. Most attempts to 
compute the downwelling longwave fluxes have relied on these retrieved profiles to furnish the 
flux using a radiative transfer model. 

Here we show results using the once daily profile contained in the data released by the 
International Satellite Cloud Climatology Project (ISCCP; Rossow and Schiffer 1991). ISCCP 
Cl data provides a daily profile of temperature and precipitable water as well as surface 
temperature at a 2.5° x 2.5° horizontal resolution. This information from TOVS (Tiros 
Operational Vertical Sounder) is used to generate clear sky downward longwave fluxes globally 
and upward fluxes only over the ocean where the once-a-day sampling is acceptable. The 
diurnal cycle of surface temperature precludes using this technique over land. The radiation 
code used to compute the clear sky downward flux is the one used in the UCLA/GLA (now 
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CSU) GCM (Harshvardhan et al. 1989; Randall et al. 1989, 1991). 

A flow diagram of the technique used to obtain surface longwave fluxes is shown in Fig. 
1. Calculations start from the relationship between the longwave CRF at the top of the 
atmosphere and at the surface (Harshvardhan et al. 1990). Results for the months of January, 
April, July and October are used in this study and represent the simulated monthly characteristics 
of this relation over the annual cycle. The longwave CRF at the top of the atmosphere as 
obtained from ERBE is combined with the relationship to obtain the longwave CRF at the 
surface for each of the four months. Then, by means of the radiative transfer model, with input 
meteorological parameters such as temperature profile and water vapor mixing ratio profile and 
sea surface temperature from ISCCP data as well as standard ozone vertical distributions, the 
clear-sky downward longwave radiative flux and the clear-sky net longwave upward radiation 
flux at the surface are obtained. Because ozone is primarily confined to the stratosphere, its 
contribution to downward longwave flux is much less than that of other parameters such as water 
vapor in the lower atmosphere. Therefore, use of a standard ozone profile is justified. 

In general, it is difficult to obtain satellite estimates of clear-sky flux that do not suffer 
from some cloud contamination due to a combination of subsensor resolution cloud elements and 
clouds that may not be detectable with the spectral intervals of the current radiometer. In this 
study, the clear sky fluxes are computed for the atmospheric structure under cloudy conditions, 
but assuming a cloud-free sky. Since the CRF at the surface and the clear-sky radiative fluxes at 
the surface are available, based on the definition of the CRF shown in equation (1), the actual 
radiative fluxes are calculated simply as the clear-sky flux plus the corresponding CRF. 

3. Results 

a. Longwave CRF at the Surface ( global ) 

The mean monthly distributions of the longwave cloud radiative forcing (CRF) at the 
surface are presented in Fig. 2 for the months of April, July and October 1985 and January 1986. 
They are derived by combining the longwave CRF at the top of the atmosphere from ERBE 
(Harrison et al. 1990) and the ratios of the CRF at the top to that at the surface. Based on the 
definition in equation (1), the longwave CRF is the difference between the mean longwave flux 
and clear sky longwave flux for the same period. Here the longwave CRFs at the top of the 
atmosphere are retrieved from ERBE satellite data according to this definition. In the original 
ERBE data, there are a few regions in the tropics where clear sky longwave radiative fluxes at 
the top of the atmosphere are unavailable because of the lack of clear sky pixels during the 
experiment period. Therefore, an interpolation is performed to make up the missing data using 
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neighboring points. 

There are several noteworthy features of these distributions. First, regions with small 
surface longwave CRF are concentrated in the tropical and the subtropical oceanic areas (Fig. 2). 
Surface CRF is small in the tropics because the boundary layer there is moist and radiatively 
opaque even for clear skies. In the central Pacific Ocean, the areas with surface longwave CRF 
below 20 Wnr 2 dominate throughout the year. This is also the case for the northern and central 
Indian Ocean in both April 1985 and January 1986. Areas with large surface longwave CRF 
occur over oceanic areas southwest of Indonesia, with a maximum of more than 60 Wm* 2 in 
October and more than 80 Wm* 2 in July 1985. Both of these maxima correspond to areas of 
tropical convective activity. In the central Atlantic Ocean, the surface longwave CRF is below 
30 Wm* 2 . 

Much larger values of surface CRFs are found over continents than over oceans during 
most of the year. The surface values in Eurasia and North American are more than 80 Wm" 2 in 
April 1985 and 60 Wm' 2 in both October 1985 and January 1986. Evidently, all these results 
correspond to persistent widespread precipitation over the continents. It is interesting to note 
that small surface CRFs (below 30 Wm" 2 ) emerge over most continents in July 1985. These low 
surface CRFs are a consequence of the large CRF ratios over continents computed from the 
GCM. The model includes the radiative effects of convective anvils but ignores shallow 
convective clouds (Harshvardhan et al. 1989). The CRF at the top is therefore quite substantial 
but the surface CRF is not. Thus the ratio is probably an overestimate. As mentioned later, a 
100 mb error in low cloud base results in a 10 Wm" 2 error in the downward longwave flux for 
complete cloud cover. If the low clouds are scattered, the error should not be significant and if 
the boundary layer is very moist, the sensitivity to cloud base is reduced further. 

There are two continents, Africa and Australia, which always have small surface CRFs 
throughout the year because of their large desert and semi-arid areas. Fig. 2 shows the surface 
CRFs there to be often below 30 Wm -2 , especially in northern Africa, where the values are 
always below 20 Wm" 2 . It is worthwhile noting that a persistent large surface CRF region exists 
in the tropical region of South America. Obviously, this results from the persistent cloud cover 
in this region. Over high latitude regions, some areas with the largest surface longwave CRFs 
are found over the region poleward of 65°S, with a maximum of 150 Wm' 2 in July and some 
areas over the Arctic, with a maximum of more than 140 Wm -2 in October and January. Two 
explanations are possible for this phenomenon. Physically, it is reasonable that persistent low 
clouds over these regions lead to high surface longwave CRF. On the other hand, an 
underestimate of high cloud over the polar regions in the GCM simulation can also result in the 
surface longwave CRFs being somewhat unreasonably high. The latter is probably the reason 
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for large longwave CRFs over high latitudes in Fig. 2. Finally, areas with the average range of 
30-50 Wm' 2 dominate most mid-latitude oceanic regions of both hemispheres. This is 
attributable to the increase in oceanic stratus in these regions. 

b. Downward Longwave Radiation Flux at the Surface ( global ) 


Downward longwave radiative fluxes at the surface for cloud- free skies, as calculated by 
the radiation code in the UCLA/GLA GCM in conjunction with the input meteorological 
parameters from the ISCCP satellite data and U. S. standard atmosphere (COESA 1976), are 
presented in Fig. 3 for April, July and October 1985 and January 1986 respectively. 

Over the oceanic areas, clear sky surface downward longwave fluxes have pronounced 
zonal distributions, especially in mid-latitudes and near polar regions. In the tropical and 
subtropical areas, the regions with large downward fluxes (larger than 400 Wm* 2 ) are centered 
over Southeast Asia in April and October 1985 and shift a little northward in July 1985 and, as 
expected, a little southward in January 1986. This is not surprising since the surface downward 
longwave fluxes for clear skies are related closely with seasonal changes in temperature and 
moisture. In the northern spring (April) and fall (October), for the quite symmetrical solar 
irradiance distribution at that time, surface downward longwave fluxes also have a symmetrical 
distribution with respect to the equator. In the northern summer (July) and winter (January), the 
maps of clear sky surface downward fluxes show a shift northward and southward respectively 
with the changing seasons. This result is consistent with the fact that clear sky downward fluxes 
are determined by the near surface temperature which is closely related to the incident solar 
energy. Evidently, the high values (larger than 400 Wm" 2 ) in Southeast Asia correspond to areas 
where, due to the high surface temperature of islands, the atmosphere is warmer than elsewhere 
along the same latitudinal belt. In particular, there is an area with values larger than 420 Wm" 2 
centered east of Papua New Guinea in January 1986, corresponding to very warm and moist near 
surface conditions. In contrast, the low values (less than 200 Wm" 2 ) near the polar regions 
correspond to areas where the atmosphere is cold and dry throughout the year. 

Over the continents, symmetrical distributions of clear sky downward fluxes disappear 
and more complicated features of distributions emerge, corresponding to the land surface. There 
is a pronounced center of low values of downward fluxes over the Tibetan Plateau in Asia 
throughout the year because the mean elevation of the Tibetan Plateau is more than 4000 m 
above sea level so that surface air temperature is quite low except in summer. However even 
though high temperatures often occur there in summer, its high elevation makes the air still very 
dry. Since the major contribution to longwave downward flux at the surface comes from the 
water vapor in the lower atmosphere, the dry air over that region results in an area of low values 
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of downward fluxes. Over southern Africa, relatively higher values along the same latitude 
result from the higher air temperature and more moist air. In northern Africa, although the 
surface temperature is very high, the dry air, due primarily to its desert areas, results in 
downward fluxes that are lower along the same latitude. In addition, a pronounced trough line 
along the Rocky Mountains in Fig. 3 represents higher altitudes and a drier atmosphere in that 
region compared to surrounding areas. 

Surface downward longwave fluxes, as derived from the formula described in equation 
(1) with given clear sky fluxes and the longwave CRF at the surface, are presented in Fig. 4 for 
April, July and October 1985 and January 1986 respectively. Over the oceanic areas, the zonal 
distribution for the clear sky case is no longer apparent. The presence of clouds increases the 
downward fluxes at the surface globally throughout the year. Fluxes with values larger than 400 
Wnr 2 are found over areas associated with the intense convective cloud systems, such as the 
ITCZ, as well as the summer and winter monsoon areas. It is reasonable because, in these 
regions, clouds occur frequently and the air temperatures are high. An area with values larger 
than 440 Wnr 2 is found over northeast Australia in January 1986, corresponding to the high 
temperatures in the austral summer. Also, an area with high values is centered over southwest 
Indonesia in July 1985. The maps for July 1985 and January 1986 compare quite well with the 
corresponding results of Wu and Cheng (1989) for 1979. Their results were obtained using a 
physical model based on HIRS 2/MSU retrievals. 

Compared with clear sky maps, larger changes occur over land than over the oceans 
throughout the year. In addition, the regions of high values are still distinguishable in these 
months except in April 1985. Low downward fluxes (less than 280 Wnr 2 ) are found over high 
latitudes and polar regions, where the precipitable water is low and air temperature is cold. As 
seen in Fig. 4, contour lines of downward fluxes are spaced very densely in high latitudes during 
the whole year. As we mentioned before, these sharp changes perhaps result from uncertainties 
in modeling the meteorological parameters in these regions, an unavoidable consequence of the 
inability to distinguish cloud cover from background snow and ice and hence compute the cloud 
forcing. In addition, seasonal variations of surface downward fluxes are still notable in Fig. 4, 
even though they are not as remarkable as in the clear sky case. 

c. Net Upward Longwave Radiation Flux at the Surface ( ocean only) 

The clear sky net upward longwave flux at the surface is the difference of the upward 
flux minus the downward surface flux assuming no clouds. The upward flux here is computed 
from the sea surface temperature, which is a reported parameter in ISCCP Cl satellite data. 
There are three surface temperatures (TS) provided by ISCCP: mean TS from a clear sky 
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composite, mean TS for IR-clear pixels, and mean TS for VIS/IR-clear pixels for only day time. 
In order to compute the surface emission, it is not necessary to use any of these fields but instead 
rely on an independent source such as the sea surface temperature provided by NOAA's Climate 
Analysis Center (CAC). This is a blend of in situ data, advanced very high resolution radiometer 
(AVHRR) satellite data, and ice data (Reynolds 1988). Fig. 5 shows the zonally averaged sea 
surface temperatures for January 1986 for the three ISCCP fields and the CAC field. In high 
latitude regions, the CAC sea surface temperature is much higher than the ISCCP values because 
ice surface temperature is set to be the freezing point of sea water (-1.8°C). Fig. 5 indicates that 
the clear sky composite TS is warmest among the three ISCCP sea surface temperatures, while 
the IR clear sky TS is the coldest. This is to be expected since IR clear pixels are contaminated 
by low level clouds. We chose the VIS/IR clear sky as the most representative quantity for this 
study. 

The difference between the surface emission computed from daytime mean TS for 
VIS/IR clear pixels and the CAC sea surface temperature is presented in Fig. 6 for January 1986. 
This difference is a measure of the uncertainty one may expect in the upward longwave radiation 
flux computed using different sources for the sea surface temperatures. The dominant feature of 
this map is a set of biases of less than 10 Wm -2 present over the tropical and subtropical regions 
except around southeast Asia, where the surface emission from mean TS for VIS/IR clear pixels 
exceeds that from CAC temperature by as much as 30 Wm" 2 . This could be a result of 
differences in the particular algorithms used by ISCCP and CAC to account for water vapor 
absorption in the atmospheric window. Also any ship observations used are not necessarily 
representative of the skin temperature of the ocean surface. The areas with differences larger 
than 10 Wm’ 2 are found over high latitudes in the Pacific and Atlantic Ocean. Positive 
differences (less than 10 Wnr 2 ) are located over the central Pacific ocean and eastern Indian 
ocean, whereas areas with negative difference are found in other locations. The maximum 
absolute difference is only around 5%. Therefore, the agreement between the two upward fluxes 
is generally good for this month. 

Fig. 7 shows the clear sky net upward longwave fluxes at the surface for April, July and 
October 1985 and January 1986 respectively. Physically the maps of clear sky net upward 
longwave fluxes primarily reflect the distribution of water vapor content in the boundary layer. 
The area with lowest values is found over Southeastern Asia throughout the year. In this region, 
even in the absence of clouds, the clear sky net upward longwave fluxes are quite small because 
of the high water vapor mixing ratio near the surface. Areas with values larger than 100 Wm* 2 
are found over the mid-latitude subsidence zones between 15° and 45° latitude in both 
hemispheres in April, October 1985 and January 1986, and in the southern hemisphere in July 
1985. These are regions in which a dry atmosphere overlies a moderately warm ocean surface. 
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Also, in July 1985, there is a minimum of less than 40 Wm' 2 off the coast of central America 
and, in the other three months, a minimum of less than 40 Wm -2 just along the west coast of 
Colombia, Ecuador and Peru; all these minima correspond to the presence of a relatively moist 
atmosphere. 

Net upward longwave flux, a difference field of the clear sky net upward flux minus the 
surface longwave CRF (Fig. 2) as described in equation (1), is presented in Fig. 8 for April, July 
and October 1985 and January 1986. In this case, the areas with smaller values are still found 
over Southeast Asia. When it is clear, the value is low because the boundary layer is moist; 
when it is cloudy, there is enhanced downward emission. Also, small net upward fluxes with 
values less than 40 Wm' 2 are found in high latitudes and polar regions where there is persistent 
cloud cover and temperatures are low throughout the year. Values larger than 80 Wm 2 are 
found over some areas in the tropical and mid-latitude zone. It is worth noting that the negative 
values near the Antarctic in July 1985 result from the large longwave CRF at the surface derived 
previously. Since negative values are possible but unlikely, one may conclude that the 
unreasonable high surface CRF is a result of the modeling and computation errors but not 
physical processes. Generally, the influence of clouds in the stratus regime causes the net 
longwave flux to be reduced by 50-60 Wm' 2 throughout the year. These maps again show 
features similar to those obtained by Wu and Cheng (1989). 

4. Discussion 

Monthly mean longwave radiation fluxes at the surface for four months have been 
determined from currently available satellite data. Because of the diurnal variation of surface 
temperature over the continents, surface net longwave fluxes, which involve the surface 
temperatures, are presented only over the oceanic areas. The method discussed here avoids the 
use of an independent estimate of the frequency of occurrence of clouds or even cloud top 
heights in determining the surface longwave fluxes. Current methods of modeling the longwave 
radiation processes in the atmosphere require certain assumptions regarding the presence of 
clouds and their horizontal and vertical extent. Because of the complicated nature of cloud 
radiation and uncertainties in the observation of clouds, large errors are inevitable in this 
procedure. The procedure used in this study provides an alternative means of obtaining 
longwave radiative fluxes at the surface without the knowledge of cloud distribution. All 
meteorological parameters required in this method can be obtained from currently available 
satellite data sets. The ISCCP data and U.S. Standard Atmosphere (COESA 1976) provide the 
profiles of temperature, water vapor mixing ratio and ozone as well as sea surface temperatures. 
The distributions of the longwave CRF at the top of the atmosphere are retrieved from ERBE 
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data which is currently being completed for several years of measurements. 

The weakest link in this procedure is the use of simulated CRF ratios that relate the CRF 
at the top and surface. Errors in the cloud generation scheme of the model used will affect the 
ratio and hence, the final product. It is also not feasible to verify these ratios observationally on 
a global scale, although it would be useful to verify the model result in some specific regions 
where simultaneous observations at the top of the atmosphere and the surface are available over 
an extended period of time. 

The final product can be compared with field data for a few regions and periods during 
which extended time observations were made from ships. Note that the fields generated are area 
averaged monthly means and comparisons with point measurements over short periods of time 
do not provide any corroboration. Observations reported by Reed and Halpem (1975) off the 
Oregon coast covered an eight week period in July - August 1973 and included measurements 
taken from two sites, one 13 km from the shore and the other 120 km away. Lind and Katsaros 
(1987) have reported measurements taken off the California coast during the first two weeks of 
November 1984. Radiation budget data from GATE (Cox and Griffith, 1979) was an amalgam 
of direct radiation measurements and modeled fields based on soundings. The authors point out 
the futility of computing area and time mean quantities directly from measurements even for a 
dedicated field campaign over a (relatively) small portion of the oceanic area of the globe. 

The daily mean net longwave radiative flux at the surface reported by Reed and Halpem 
(1975) varied from 71 Wm' 2 for days when the daily mean cloud cover was 10% - 20% to 
between 11-15 Wm -2 for 70% - 100% cloud cover. An unweighted mean of the ten daily mean 
values taken at two stations over the period July 5 - August 26 is 33 Wm' 2 . The corresponding 
estimate for July 1985 from Fig. 8 is around 40 Wm' 2 . This suggests that the procedure used 
here is able to provide a good measure of the mean cloud fraction which is the strongest 
determinant of the net longwave radiation for this region. This bears out our thesis that a reliable 
measure of the effect of clouds is a necessary condition for obtaining global fields of the net 
longwave radiation at the surface. 

Lind and Katsaros (1987) have reported ship based observations of the upward and 
downward longwave radiation taken off the California coast from October 30 - November 14 
1984. The daily mean net longwave radiation from R/P FLIP ranged from 1 1 Wm' 2 to 69 Wm' 2 . 
Although there is no reported cloud cover, inspection of the daily mean insolation shows that the 
extremes correspond to days of complete cloud cover and essentially clear skies, respectively. 
The mean for the 15 day period at the single station was 45 Wm* 2 . The mean for October 1985 
and January 1986 from Fig. 8 is 75 Wm' 2 - This compares favorably with the satellite and model 
based results of Wu and Cheng (1989) for January 1979. The discrepancy with the point 
measurement is significant but it is difficult to draw any conclusions based on the information 
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available. From Fig. 3, the clear sky downward longwave flux for this region is about 315 Wm' 2 
in October 1985 and 275 Wnr 2 in January 1986. The data from Lind and Katsaros (1987) 
indicate that the minimum daily mean value over the period was 340 Wm* 2 , presumably on the 
clearest day. Inspection of the temporal data shows that the absolute minimum value reached 
was 300 Wm" 2 on the night of November 7-8. This indicates that the clear sky element of the 
procedure is acceptable. The discrepancy is in the cloud estimate and there is no reasonable 
means of comparing the point measurement with an area average (for a different year). This 
highlights the difficulties inherent in validating global fields for a quantity that is measured only 
occasionally. 

Some generalizations can be made about the uncertainty in the global fields. Fung et al. 
(1984) state that a 1 g kg" 1 uncertainty in the water vapor mixing ratio at low levels and a 100 
mb error in the cloud base each translate to a 10 Wm" 2 uncertainty in the net longwave flux over 
oceans. The presence of low level haze which cannot be detected from space and is not 
considered in the GCM is equivalent to underestimating low cloud cover, hence surface forcing. 
In certain regions this error could range from 5-10 Wm' 2 . Moreover, if one makes the gross 
assumption that low clouds cover the oceans 50% of the time everywhere, it may be shown that 
the net longwave at the surface will range from 40 - 60 Wm' 2 for any reasonable temperature 
and mixing ratio profile. The departure from this range of values shown in Fig. 8 is the true 
information contained in the maps and reflects the satellite based cloud information that has been 
used. 

5. Conclusions 

An attempt has been made to produce monthly mean global fields of the net longwave 
radiation flux at the surface over the oceans without resorting to direct measurements. The key 
ingredients in the technique are the satellite derived temperature and moisture profiles (which are 
available operationally), the top of the atmosphere cloud radiative forcing from ERBE (which is 
an experimental data product but could be available in the future) and the cloud distributions 
from a general circulation model. It should be stressed that the actual cloud cover generated by 
the model (which is subject to a great deal of uncertainty) is not used directly, but only 
information on cloud type is used through a ratio of the cloud radiative forcing at the top and 
surface. This parameter has the advantage that it involves simulated radiation fluxes and not the 
cloud fraction. The latter quantity may vary considerably from model to model and also be quite 
different from satellite derived estimates. However, the radiation parameterization is usually 
such that the fluxes at the top of the atmosphere simulated by these models compare quite 
favorably with observations. Examples of this apparent contradiction are in Harshvardhan et al. 
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(1989), and Kiehl and Ramanathan (1990). 

Although it is preferable to map global fields of the longwave radiation at the surface 
using space and ground based measurements alone, it is evident that the introduction of large 
scale numerical models into this effort is unavoidable. The surface is inaccessible to space based 
instruments at these wavelengths under cloudy conditions. Moreover, since the fluxes are 
extremely sensitive to water vapor mixing ratios near the surface, even clear sky estimates are 
subject to a great deal of uncertainty. Any global surface measurement program over the oceans 
is impractical. Hybrid techniques such as the one reported here that use several sources of data, 
both real and simulated, are the only options. 
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Figure Captions 

Figure 1. Flow diagram of the procedure to obtain maps of the monthly mean net upward 
longwave radiation flux at the surface using information provided by a general 
circulation model (GCM), data from the Earth Radiation Budget Experiment 
(ERBE) and the International Satellite Cloud Climatology Project (ISCCP). 

Figure 2. Monthly mean longwave cloud radiative forcing at the surface for April, July, 
October 1985 and January 1986 obtained from ERBE top-of-the-atmosphere 
cloud forcing and GCM simulations of cloudiness. 

Figure 3. Monthly mean clear sky downward longwave radiation flux at the surface for 

April, July, October 1985 and January 1986 obtained from TOVS profiles on the 
ISCCP Cl tapes and a broad band radiation code. 

Figure 4. Monthly mean atmospheric downward longwave radiation flux at the surface for 
April, July, October 1985 and January 1986 obtained from the clear sky values 
shown in Figure 3 and the cloud forcing shown in Figure 2. 

Figure 5. Zonal mean sea surface temperatures (SST) for January 1986 from three different 
parameters on the ISCCP Cl tapes and the blended SST from the Climate 
Analysis Center (CAC). 

Figure 6. Difference between the monthly mean surface emission in Wm' 2 for January 
1986 computed using the CAC SST and the VIS/IR clear TS on the ISCCP Cl 
tapes. A positive difference indicates that the surface emission implied by the 
CAC SST is higher. 

Figure 7. Monthly mean clear sky net upward longwave radiation flux at the ocean surface 
for April, July, October 1985 and January 1986. 

Figure 8. Monthly mean atmospheric net upward longwave radiation flux at the ocean 
surface for April, July, October 1985 and January 1986. 
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Clear Sky Net Upward Longwave Flux at the Surface 
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ABSTRACT 


Computational results have been obtained for the spherical albedo, global transmission, 
and global absorption of plane-parallel layers composed of cloud droplets. These computations, 
obtained using the doubling method for the entire range of single scattering albedos (0 < a> 0 < 1) 

and for optical depths between 0.1 and 100, are compared with corresponding results obtained 
using selected multiple scattering approximations. Both the relative and absolute accuracies of 
asymptotic theory for thick layers, three diffuse two-stream approximations, and two integrated 
two-stream approximations are presented as a function of optical thickness and single scattering 
albedo for a scattering phase function representative of cloud droplets at visible wavelengths. 
The spherical albedo and global absorption computed using asymptotic theory are found to be 
accurate to better than 5% for all values of the single scattering albedo, provided the optical 
thickness exceeds about 2 . The diffuse two-stream approximations have relative accuracies that 
are much worse than 5% for the spherical albedo over most of the parameter space, yet are 
accurate to within 5% in the global absorption when the absorption is significant. The integrated 
delta-Eddington scheme appears to be the most suitable model over the entire range of variables, 
generally producing relative errors of less than 5% in both the spherical albedo and global 
absorption. 
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Introduction 


The role of clouds in determining the Earth’s radiation budget has led to increased 
interest in the parameterization of the radiative properties of cloud layers in numerical 
atmospheric models. Recent work has been concerned with relating cloud microphysics to 
optical properties (Slingo 1989) which can then be used in radiative transfer schemes within 
models. Most models now use some form of approximation to compute cloud radiative 
properties, such as the plane albedo from a given set of optical properties (optical thickness, 
single scattering albedo, etc.) Whereas in the past these optical properties were generally fixed, 
there is now increasing use of interactive schemes in which cloud optical properties are 
generated internally by the model (Charlock and Ramanathan 1985; Harshvardhan et al. 1989). 

As cloud fields evolve during a model integration, the optical properties of the generated 
clouds and models of gaseous absorption are used in a radiative transfer scheme to provide the 
shortwave and longwave radiative energy field through the atmosphere. These computations 
need to be carried out at each model grid point at least every time the model cloud fields are 
updated. In models that resolve the diurnal cycle, this could be every three hours of simulated 
time or even hourly. The computational burden is such that rapid yet accurate techniques are 
essential. In the shortwave, a common procedure is the computation of cloud layer properties by 
a two- stream method and the adding of radiative fluxes through the atmosphere in an energy 
conserving scheme (Lads and Hansen 1974; Coakley et al. 1983; Charlock and Ramanathan 
1985; Geleyn and Hollingsworth 1979; Harshvardhan et al. 1987) although the two-stream 
equations can also be solved directly for multiple layers using matrix solvers (Wiscombe 1977; 
Toon et al. 1989). The flux adding method is essentially a severely truncated form of the 
adding-doubling method (Hansen and Travis 1974) using upward and downward fluxes instead 
of intensities. 

In order to compute radiative fluxes through several atmospheric layers by the flux 
adding method, the radiative properties of cloud layers for two different sources are required 
(Harshvardhan et al. 1987; Kiehl et al. 1987). When collimated solar radiation is incident on an 
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isolated cloud layer at some zenith angle with respect to the vertical direction, the fluxes 

emergent from the layer in the upward and downward directions are determined by the plane 

albedo and total transmission of the layer respectively. If the incident source is diffuse, the 

emergent flux may be obtained by an angular integration over the incident intensity field. In two- 

stream methods, the angular distribution of the incident intensity field is not resolved and a 

common practice is to assume an isotropic diffuse source. For example, in a multi- layer cloud 

system, the diffuse solar flux transmitted through the upper layer is the incident source for the 

lower layer. Also in the case of a cloud layer overlying a reflecting ground surface, multiple 

reflections between the cloud and ground are considered by assuming an isotropic diffuse source 

at the bottom boundary of the cloud layer. These diffuse radiative properties have also been 

/ 

used in the past to provide estimates of global effects of aerosol layers (Chylek and Coakley 
1974). A comprehensive study of the accuracy of various multiple scattering approximations for 
the plane albedo, total transmission and fractional absorption of isolated cloud layers 
corresponding to incident collimated radiation was presented by King and Harshvardhan (1986a, 
b). The present study complements the earlier one in assessing the accuracy of various 
approximations for calculating the radiative properties of cloud and aerosol layers for an incident 
isotropic diffuse source. 

The presentation follows the organization of King and Harshvardhan (1986a; hereafter 
referred to as KH). Section 2 discusses multiple scattering calculations used to obtain the diffuse 
radiative properties of cloud layers of varying optical thicknesses and single scattering albedos. 
These computational results, obtained with the doubling method, will be considered the 
benchmark solutions with which various multiple scattering approximations will be compared. 
Section 3 introduces the asymptotic theory approximation and the general class of two-stream 
approximations that we will consider. Section 4 presents the results of the comparison between 
the approximate and exact results in terms of absolute and relative differences. A discussion of 
the results follows in Section 5. Section 6 is a summary including recommendations for using 
these approximations. 
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2. Multiple scattering computations 

To assess the accuracy of various multiple scattering approximations, radiative transfer 
computations were performed using the doubling method described by Hansen and Travis 
(1974), together with the invariant imbedding initialization described by King (1983). These 
computations were performed for a cloud drop size distribution typical of fair weather cumulus 
(FWC) clouds (Hansen 1971), and were performed at a wavelength X = 0.7 54pm assuming a 
refractive index of liquid water m = 1.332. A detailed description of the cloud model, together 
with an illustration of the single scattering phase function, can be found in KH. The azimuth- 
independent terms of the reflection and transmission functions were used to obtain the plane 
albedo r(x t ,p 0 ) and total transmission t(x t ,p 0 ) as a function of x„ the total optical thickness of the 
layer, and p 0 , the cosine of the solar zenith angle. In terms of these functions the spherical 
albedo, global transmission and global absorption of the layer are given by 

1 

r(x t ) = 2jr(x li p 0 )p 0 dp 0 , 

0 

l 

t(x t ) = 2 Jt(x t> p 0 )p 0 dp 0j 
0 

a(x) = 1 - r(x t ) - t(x t ). 

In order to cover a wide range of applications, these computations were performed for 
values of the single scattering albedo ranging from pure absorption (co 0 = 0) to conservative 
scattering (to 0 = 1). The single scattering phase function was left unchanged such that all 

computations apply to a phase function having an asymmetry factor g = 0.843. 

Figure 1 illustrates numerical computations of the spherical albedo [f(x t )], global 
transmission [t(x t )] and global absorption [a(x t )] as a function of co 0 and x t . The doubling 


( 1 ) 

( 2 ) 

(3) 
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computations used to generate these results were obtained at twelve optical depths 0.0625, 0.125, 
..., 128 interleaved with another set of eleven optical depths 0.0884, 0.1768, ..., 90.51. Each set 
of doubling computations was itself made at each of 31 values of the single scattering albedo. 
The single scattering albedo scale is linear in the similarity parameter s, defined by 


s 


1 - COn 


1 - “og. 


1/2 


(4) 


This makes it possible to expand the scale in the vicinity of conservative scattering (co 0 = 1) and 
still to span the full range 0 < to 0 < 1. The angular computations, including the integration in (1) 
and (2), were performed at 80 Gaussian quadrature points. As in KH, the computed results were 
first interpolated to generate a 300 x 300 matrix prior to plotting. The interpolated arrays 
represent the exact results to which the radiative transfer approximations are compared in 
Section 4. 

It is perhaps pertinent to point out certain features of the radiative properties illustrated 
in Figure 1. For conservative or very weakly absorbing layers, the spherical albedo increases 
rapidly with increasing optical thickness for small values of x t and then much more slowly as T t 

becomes large. This is the well known non- linear behavior that leads to problems in estimating 
area- averaged albedos for a non- homogeneous cloud layer (Harshvardhan and Randall 1985). 
For moderate to strong absorption, the saturation of both the spherical albedo and global 
absorption at optical thicknesses of about 10 or even less is the most striking feature of Figure 1. 
In the near infrared, this implies that cloud absorption is primarily a function of the single 
scattering albedo and not the optical thickness once the cloud layer is several hundred meters 
thick (Twomey 1976). The importance of determining the spectral dependence of d) 0 for cloud 

layers and the development of accurate parameterizations for inclusion in radiative transfer 
models follows from this observation (King et al. 1990; Fouquart et al. 1991). 
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3. Radiative transfer approximations 

Three classes of approximations will be considered here for comparison with the multiple 
scattering results presented above. In all cases, analytic or easily integrable functions relate the 
radiative properties to the optical properties. The three approximations we will consider are 
asymptotic theory for thick layers, diffuse two-stream approximations, and integrated two- 
stream approximations. Although there are several variations of two-stream approximations, 
only a few common and representative models will be considered. 
a. Asymptotic theory 

Asymptotic theory is a rigorous solution to the equation of transfer in optically thick 
layers, and as such, makes no assumption about the angular distribution of scattered radiation 
within the medium. Expressions for the plane albedo and total transmission of an optically thick 
layer under collimated illumination conditions can be found in KH and will not be repeated here. 
From these expressions it can be shown that the asymptotic theory approximations for the 

A A A 

spherical albedo [r(x t )], global transmission [t(t t )] and global absorption [a(x t )| are given by 


r(x t ) = 

mn 2 le' 21cTt 
r °° — „ -?W ’ 

1 - l 2 e 2kTl 

(5) 

n 

/*> 
4— » 

mn 2 e" kTt 

(6) 

1 - l 2 e ' 2kTt ’ 

A 

a(x t ) = 

A A 

1- r(x t )- t(T t ) , 

(7) 


for nonconservative scattering (co 0 < 1). In these expressions r^ is the spherical albedo of a 

semi-infinite atmosphere and m, n, 1 and k are constants (coefficients) that depend primarily on 
the similarity parameter given by (4). All of the functions and constants that appear in these 
expressions can be computed by equating asymptotic formulae and doubling results at three 
values of the optical thickness for which asymptotic theory is valid (viz., x t = 8, 16 and 32), as 
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first pointed out by van de Hulst (1968). Similarity relations for calculating r^ (denoted A* by 

van de Hulst 1968), m, n, 1 and k as a function of s for the full range 0 < s < 1 can be found in 
Table 1 of King et al. (1990) Once these coefficients have been computed, expressions for all 

of the radiative properties are analytic functions that can be computed rapidly within a radiative 
transfer code. 

For the special case of conservative scattering (o) 0 = 1), Eqs. (5) and (6) reduce to 


f (x t ) = 1 - 


4 

3(1 - g)(x t + 2q 0 ) ’ 


( 8 ) 


6 4 

t(X ' ) '3(l-g)(T t + 2q 0 ) 


(9) 


where q 0 is the extrapolation length. The reduced extrapolation length q' = (l-g)q 0 is known to 
range between 0.709 and 0.715 for all possible phase functions (van de Hulst 1980), and has the 
value q' = 0.715 for the phase function used here. Again, one is left with simple analytic 

A A 

functions describing the variation of r(x t ) and t(x t ) as a function of x t for a given asymmetry 
factor g. The set of equations (5) - (9) form the approximations for the diffuse radiative 
properties of a medium based on asymptotic theory. 
b. Diffuse two-stream approximations 

In the absence of any direct collimated beam, the two- stream equations of radiative 
transfer result in a set of differential equations for the upward and downward diffuse fluxes 
'F^(x) (Coakley and Chylek 1975; Meador and Weaver 1980): 


d F-QO 

dx 


= 7 1 F'(x)- Y 2 f+ ( x )’ 


dF+(x) 


= Y 2 f '(x) - y,F + (x), 


( 10 ) 


dx 


( 11 ) 
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where F'(x) represents the upward flux and F + (x) the downward flux at optical depth x. The 
equations can easily be solved subject to the boundary conditions 


f+(0) = f 0 , 

(12) 

o' 

ii 

4-» 

Uh 

(13) 


for a diffuse isotropic source incident at the top boundary of the layer (or cloud) and for which 
no illumination is incident from below. The spherical albedo is thus obtained from the 
expression 


f(t t )=F-(0)/F 0 


and the global transmission from 


t(T t )=F+(T t )/F 0 , 


(14) 

(15) 


For nonconservative scattering (co 0 < 1), the solution may be obtained in the form (Coakley and 
Chylek 1975; Meador and Weaver 1980 ) 



■fed - e' 2 * 1 ) 

k + Y, + (k-y,) e" 2ktt 


too = 


2k e' kTl 


k + Yi + (k-Yj)e 


■2kt t ’ 


(16) 


(17) 


and for conservative scattering (co 0 = 1): 



7i x t 

1+YiTt ’ 


(IS) 
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A A 

t (T t ) = 1 - r (x t ) . (19) 

In (16) - (19), the coefficients Yl and y 2 depend on the particular two-stream approximation, with 
the diffusion exponent k defined as 


k = ( Yl 2 


- y 2 2 ) ,/2 - 


( 20 ) 


Table 1 lists three diffuse two-stream models used for this study and the corresponding 
values of Yl and y 2 . The discrete ordinates model is identified as the quadrature scheme by 

Meador and Weaver (1980) and Toon et al. (1989). The hemispheric mean model defined by 

Toon et al. (1989) is similar to the Coakley-Chylek model II referred to by KH and first 

introduced by Chylek and Coakley (1974). The two stream model used by Sagan and Pollack 

(1967) has coefficients similar to those of both of the above mentioned models. Instead of the 

asymmetry parameter g, some two stream models use the average backscatter fraction (3, which 

is defined in KH and readily computed from the backscatter fraction (3(p. 0 ), introduced by 

Coakley and Chylek (1975) and Zdunkowski et al. (1980) to compute the radiative properties of 

layers for collimated incident sources. The Eddington model has, of course, been used widely 

(Shettle and Weinman, 1970). The set of equations (16) - (20) are used to compute the diffuse 

radiative properties for the two-stream approximations. It should be noted that these expressions 

have fairly simple analytic forms that favor rapid computation. 

c. Integrated two-stream approximations 

Extensive discussion of two-stream approximations for a collimated source can be found 

in KH as well as in earlier work, in particular the comprehensive treatment by Meador and 

a . . 

Weaver (1980). Expressions for the approximate plane albedo [r(x t , |i 0 )], total transmission 

A A 

[t(x t , (L 0 )] and fractional absorption [a(x t ,|i 0 )] are the set of Eqs. (21) - (29) in KH. These 
expressions include the transformations that are required in the case of delta-scaling (Joseph et 


9 


al. 1976). To obtain comparable expressions for the diffuse radiative properties, r(t t ,|i 0 ) and 

A 

t(t t , fi 0 ) must be integrated according to Eqs. (1) and (2). However, these expressions are quite 
complicated and thus integration in a closed form is not generally practical. 

An analytic expression for the spherical albedo in the Eddington and delta-Eddington 
approximations has been obtained by Wiscombe and Warren (1980) and involves exponential 

A 

integrals that are not conducive to rapid computation within a model. For this study, r(T t ,Ho) and 

A 

t (t,,p 0 ) obtained by the delta-Eddington approximation were numerically integrated to provide 

A A 

r(x t ) and t(t t ). KH found that the delta-Eddington approximation for collimated illumination 
conditions is quite accurate over a wide range of and |j. 0 especially when cd 0 is near unity. A 
model that performs well for optically thin layers over the limited range of oo 0 studied by KH is 

/ f 

the plane albedo scheme of Coakley and Chylek (1975), designated Coakley-Chylek model I by 
KH. Two-stream methods for collimated sources require a third coefficient, y 3 , which appears 
with the source term and is thus not included in Eqs. (10) and (11). The expressions for y 3 used 
by the two integrated models presented here are given in Table 1. 

The integrations in Eqs. (1) and (2) required to obtain the diffuse properties are 
performed using 80 point Gaussian quadrature and the results should be considered identical to 
an analytic solution for all practical purposes. The general form of the quadrature summation is 

N 

r (T t ) = 2 X r Cq.^PiW; (21) 

i=l 

where (ij are the Gaussian quadrature points on the half space and w ; are the corresponding 
Gaussian weights. However, this detailed integration is of no practical value because the 
computational burden is onerous when applied to a global climate model. We have therefore 
also included results for the delta-Eddington and Coakley-Chylek (I) models integrated using 
two-point and four-point quadrature, respectively. The diffuse radiative properties can then be 
obtained with a computational effort comparable to that required to compute properties for 
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collimated radiation. 

4. Results 

We have examined both the absolute and relative accuracies of the spherical albedo, 
global transmission and global absorption as a function of x t and co 0 for the asymptotic 

approximation as well as the Eddington, discrete ordinates and hemispheric mean diffuse two- 
stream approximations. Other diffuse two-stream approximations that we have examined 
generally yield somewhat poorer results when compared to our doubling benchmark 
calculations. In addition, we have considered the integrated delta-Eddington and 
Coakley-Chylek (I) approximations computed using both 80 points and a limited number of 
Gaussian quadrature points for integration over the solar zenith angle. 

Figure 2 illustrates a 4 x 3 plot composite of results for the absolute difference in the 
spherical albedo, global transmission and global absorption for four of these models, where the 
first row applies to asymptotic theory and succeeding rows to the Eddington, discrete ordinates 
and hemispheric mean approximations. Individual plots in the first column of Fig. 1 represent 
absolute errors in the spherical albedo, defined as 

A 

Ar(t t , co 0 ) = r (x t , co 0 ) - r(x^ co 0 ), (22) 

with succeeding columns representing corresponding errors in global transmission [At(x t , o) 0 )] 
and global absorption [Aa(x t , co 0 )]. With these definitions, positive (negative) errors indicate that 
the radiative transfer approximation overestimates (underestimates) the exact solution, taken as 
the computational results presented in Fig. 1. The relative errors in the spherical albedo, global 
transmission and global absorption are presented in Fig. 3, and are given in percent. It is 
necessary to consider the performance of a particular model in both a relative and an absolute 
sense in order to delineate a range of acceptability. 

Individual contour plots in Figs. 2 and 3 have been shaded to draw attention to those 
regions of greatest accuracy. For example, asymptotic theory is seen to be accurate to within 5 % 
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in reflection and absorption for x t > 2 and for all values of oo 0 . In transmission, relative errors 

exceed 5 % for co 0 < 0.90 and 2 < x t < 8 but the absolute errors are so small (<0.03) that the 

approximation could probably still be used without serious adverse results. It is evident from 

Figs. 2 and 3 that the asymptotic approximation provides accurate results for all three diffuse 
radiative properties over the entire range of co 0 as long as x t > 2. 

The three diffuse two-stream models considered here are seen to yield unacceptable 
errors in one or more of the radiative properties over regions that would normally be 
encountered in modeling applications. Although the range of acceptability will depend on the 
particular application, one can consider a 5% error in the spherical albedo as a standard for 
comparison. The spherical albedo is usually the parameter of choice in estimating the sensitivity 
of any radiative perturbation. However, when the value itself is small, an absolute error criterion 
is more useful. For optically thin layers, the absolute errors in spherical albedo are generally 
less than 0.01 for the discrete ordinates and hemispheric mean approximations. Errors in global 

transmission are similar for all three models while the Eddington and hemispheric mean models 
are successful in estimating the global absorption of a layer when co 0 > 0.99 and x t < 10 with 

errors of less than 1%. If the range of acceptability is relaxed to 5 %, then the Eddington and 
hemispheric mean models can be used for absorption when cd 0 is as low as 0.95 except for 

optically thick layers. This covers the range of single scattering albedo encountered in water 
clouds throughout the visible and near-infrared spectrum (King et al. 1990). 

The two integrated two-stream methods studied in this investigaton provide more 
accurate results for all three diffuse radiative properties as shown in Figs. 4 and 5. The delta- 
Eddington model was shown by KH to be highly successful in estimating the plane albedo for 
conservative scattering. There was a marked degradation of performance when nonconservative 
cases were considered. The present study shows that this model, when integrated over an 
isotropic diffuse incident source, provides excellent results for the spherical albedo and global 
transmission over most of the range of t t and co 0 . Errors in excess of 10% in global absorption 
are present for moderate optical depths (0.5 < x t < 5) when cd 0 exceeds about 0.95. However, it 
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may be seen from Fig. 4 that the absolute errors in global absorption are less than 0.02 

throughout this region. In addition, the large relative error in global transmission for optically 

thick absorbing layers is irrelevant since the global transmission is itself close to zero as is the 

/ 

absolute error. The Coakley-Chylek (I) model provides results of comparable accuracy for 
optically thin layers. This is not surprising since KH showed that it was the most accurate of the 
two-stream models for this case. However, the delta-Eddington model when integrated over all 
incident angles is nearly as accurate as the Coakley-Chylek (I) model for optically thin layers. 
Moreover, the accuracy of the integrated delta-Eddington model does not degrade as rapidly at 
higher optical depths. 

As mentioned previously these two models would only be of academic interest if a 
rigorous numerical integration was required for every computation of the diffuse radiative 
properties. We have therefore also presented results obtained using a limited number of 
quadrature points in the integration over solar zenith angle [cf. Eq. (21)]. As can be seen from 
the second panel of Figs. 4 and 5, a two-point integration of the delta-Eddington models yields 
accuracies that are comparable to the accuracy obtained using an 80 point integration. However, 

f 

for the Coakley-Chylek (I) model it is necessary to use a four point integration to obtain results 
that are of comparable accuracy. 

5. Discussion 

Although the results presented here are not exhaustive in the sense that all possible 
approximations -have not been tested, we feel they are representative of what one 'might expect 
for any class of model. All the schemes are computationally efficient and it is not necessary to 
perform a rigorous integration for the models based on incident collimated sources. The 
approximations presented here can be incorporated into a multi-layer radiative transfer module 
that may be added to the radiation code in a numerical model. 

All computations presented here were obtained for a FWC drop size distribution having 
an asymmetry factor g = 0.843. Variations along the co 0 axis can therefore be viewed as 

representing the effect of altering the gaseous absorption in the layer at a particular wavelength 
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or to some extent variations in the wavelength if g does not vary too greatly. This would cover 
the solar near-infrared spectrum over which (1 - co 0 ) varies by several orders of magnitude while 

g generally lies between 0.80 and 0.90 (cf. King et al. 1990). 

As found by KH for collimated radiative properties, the asymptotic approximation yields 
consistently excellent results for optically thick layers, regardless of single scattering albedo and 

solar zenith angle. Figures 2 and 3 show that the same is true for the diffuse radiative properties 
as long as x t >2. In a numerical model with internally generated cloud optical properties, this 

requirement will not always be met. Errors become unacceptably large when t t < 1. For this 
reason, the asymptotic approximation should only be used when it is known a priori that x t > 2 at 

all times. This is the one serious shortcoming of an otherwise simple and accurate model. The 
method also requires a pre-computed table of coefficients m, n, k, 1 and r^ or analytic forms that 

compute these quantities within the program. Analytic expressions for these coefficients in 
terms of the similarity parameter can be found in King et al. (1990), which further discusses a 
remote sensing application of asymptotic theory. 

The three diffuse two-stream models presented here are the simplest to implement in a 
numerical atmospheric model and are the most computationally efficient, but their accuracy is 
limited to certain regions of the parameter space. They are also not uniformly accurate for all 
three radiative properties. This is especially true in the Eddington approximation, where the 
spherical albedo is frequently too inaccurate to be of any value in a numerical model. In 
addition, the Eddington model yields unphysical values of the spherical albedo and global 
absorption when absorption is very large (King and Harshvardhan 1986b). This situation arises 
occasionally in the water vapor bands of the near-infrared and frequently in the thermal infrared. 
The problem can be rectified in a computer code with the addition of a check for unphysical 
values which could then be forced to the condition of zero reflection. The discrete ordinates 
model does not suffer from this limitation and generally yields better results for the spherical 
albedo than does the Eddington approximation. The somewhat poorer results for global 
absorption are not too important since the absolute errors are small in this case. The hemispheric 
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mean model yields results very similar to the discrete ordinates model except for global 
absorption. The smaller relative errors for weak absorption is an especially attractive feature of 
the hemispheric mean model which otherwise suffers from the fact that it tends to overestimate 

i 

the spherical albedo by more than 5 % for the very important case of nearly conservative 
optically thick layers. 

The integrated delta-Eddington model yields excellent results for all three radiative 
properties over the entire range of optical properties that are encountered in the radiation code of 
a numerical atmospheric model. In fact, errors in the diffuse radiative properties are generally 
smaller than the errors found by KH for collimated radiative properties, with no unphysical 
results anywhere in the parameter space. There has obviously been some cancellation of errors 
in the angular integration. As mentioned earlier, the one error-prone region is moderate optical 
thickness and weak absorption. This was also true for the errors in fractional absorption for a 
collimated source. Since the direct beam is usually handled by a delta-Eddington or similar 
approximation, the coefficients and functions used for this model are usually already present in a 
numerical model. There is however an extra computational overhead in the angular integration 
in that planar properties need to be computed at several angles and then numerically integrated. 
However, as seen in Figs. 4 and 5, these computations need be carried out at only two points to 
yield results comparable to a detailed numerical integration. 

The integrated Coakley-Chylek (I) model is of limited value, except perhaps for optically 
thin weakly absorbing layers. There is also an added computational burden since at least four 
angular computations are required for the phase function used here. For collimated radiative 
properties and for optically thin layers, KH found that this model was superior to the delta- 
Eddington model. For diffuse radiative properties, on the other hand, we find that there is little 
advantage in using the Coakley-Chylek (I) model, even for optically thin layers. 
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6. Summary and Recommendations 

In the present study the spherical albedo, global transmission and global absorption 
computed by various radiative transfer approximations have been compared with doubling 
computations as a function of optical thickness and single scattering albedo. Since the entire 
range of co 0 has been considered for optical depths from 0.1 to 100, the results presented here 

can be utilized to decide which approximate method is the most accurate for a particular 
application. The results presented here should be considered in parallel with the findings of KH 
regarding the plane albedo, total transmission and fractional absorption for a collimated incident 
source. 

In order to summarize the results of this study, it is useful to present composite figures 
extracted from the individual figures to highlight regions of highest accuracy. Following van de 
Hulst (1980), we show in Fig. 6 the regions for which a particular model is accurate to within 
1 % and 5%. Only those models that are reasonably accurate in the particular radiative property 
have been included. These models include asymptotic theory, the two-point delta-Eddington 
method and the four-point Coakley-Chylek (I) method. Although the hemispheric mean model 
yields acceptable results for the global absorption, it is not included here because results for the 
spherical albedo are generally poor. 

At the 1% ( 5 %) level, asymptotic theory can be used for all co 0 as long as 
T t > 3.5 (2). For smaller optical depths, there is a choice that can be made between the delta- 

Eddington and Coakley-Chylek(I) models, but our recommendation is to use the delta-Eddington 
method. Many general circulation models are already using this method to compute collimated 
radiative properties and the additional overhead incurred in the two-point integration should be 

minimal. If a scheme is needed to span the entire domain, the asymptotic method should not be 
used since its performance deteriorates very rapidly for x t < 3. For this situation, typical of 

GCM applications, the integrated delta-Eddington scheme should yield acceptable results. 

The overall errors for a multi-layer cloud system over a reflecting surface will depend on 
the optical thickness and single scattering albedo of the individual layers. At present it is felt 
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that errors in parameterizing the band averaged single scattering albedo of cloud layers in the 
near-infrared will dominate errors in approximating the radiative properties of individual layers 
(Fouquart et al. 1991). For example, the use of a single value of ft) 0 to represent the entire solar 

near infra-red can result in errors in the layer absorption of several hundred percent (Slingo 
1989). The sensitivity of all radiative properties to co 0 can be appreciated by inspection of 

Figure 1. Since any scheme has to limit the number of bands for computational efficiency, the 
selection of these bands and the average absorbing properties used could determine the overall 
accuracy. However, for a given set of x t and co 0 the results presented in this study could act as a 

guide for choosing an appropriate model. Finally, it is pertinent to mention that these accuracies 
refer to an idealized plane parallel model. There is, of course, the additional problem of 
representing inhomogeneous cloud systems including geometric effects (Harshvardhan and 
Thomas 1984; Stephens 1988), a problem not considered in this study. 
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Table 1. Summary of yj coefficients in selected two-stream approximations. 


Method 

Yi 

Y 2 

y 3 

Diffuse 

Eddington 

1/4 [7- (0 0 (4+3g)] 

- 1/4 [1- G3 0 (4— 3g)] 

— 

discrete ordinates 

\f3/2[2- co 0 (l+ g)] 

V3/2[co 0 (l-g)] 

— 

hemispheric mean 

2- co 0 (l+g) 

®o(l- g) 

— 

Integrated 

delta- Eddington 

1/4 [7- cOo(4+3g / )] 

- 1/4 [1- C0o(4-3g')] 

1/4(2- 3gV 0 ) 

Coakley-Chylek (I) 

{1- co 0 [l- p(M- 0 )]}/ki 0 

CDoPOloVUo 

P(Ho) 


co 0 ' = (1- g 2 )oy(l- co 0 g 2 ) 
g' = g/(l+g) 
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FIGURE CAPTIONS 

Fig. 1. Doubling computations of the (a) spherical albedo r(x t , co 0 ), (b) global transmission t 
(x t , co 0 ) and (c) global absorption a(x t , co 0 ) as a function of optical thickness and 
single scattering albedo for the FWC phase function. The single scattering albedo 
scale is linear in the similarity parameter, defined by Eq. (4). 

Fig. 2. Absolute accuracy of asymptotic theory, Eddington, discrete ordinates and 
hemispheric mean approximations to the spherical albedo, global transmission and 
global absorption as a function of optical thickness and single scattering albedo. The 
FWC phase function is assumed throughout. 

Fig. 3. As in Fig. 2 except for relative accuracies (in percent). 

Fig. 4. As in Fig. 2 but for the integrated delta-Eddington, two-point integrated delta- 
Eddington, integrated Coakley-Chylek (I) and four-point integrated Coakley- Chylek 
(I) approximations. 

Fig. 5. As in Fig. 4 except for relative accuracies (in percent). 

Fig. 6. The domain of validity of selected approximations to the diffuse radiative properties 
of cloud layers as a function of optical thickness (x t ) and single scattering albedo 
(to 0 ). The upper panels correspond to a relative accuracy of 1% and the lower panels 
to 5%. The single scattering albedo scale is linear in the similarity parameter. The 
FWC phase function is assumed throughout. 
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